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Abstract — In this paper, we investigate synchronization in 
a small-world network of coupled nonlinear oscillators. This 
network is constructed by introducing random shortcuts in a 
nearest-neighbors ring. The local stability of the synchronous 
state is closely related with the support of the eigenvalue distri- 
bution of the Laplacian matrix of the network. We introduce, 
for the first time, analytical expressions for the first three 
moments of the eigenvalue distribution of the Laplacian matrix 
as a function of the probability of shortcuts and the connectivity 
of the underlying nearest-neighbor coupled ring. We apply these 
expressions to estimate the spectral support of the Laplacian 
matrix in order to predict synchronization in small-world 
networks. We verify the efficiency of our predictions with 
numerical simulations. 

I. Introduction 

In recent years, systems of dynamical nodes intercon- 
nected through a complex network have attracted a good deal 
of attention [20]. Biological and chemical networks, neural 
networks, social and economic networks [9], the power grid, 
the Internet and the World Wide Web [8] are examples of 
the wide range of applications that motivate this interest 
(see also [15], [4] and references therein). Several modeling 
approaches can be found in the literature [8], [22], [1]. In 
this paper, we focus our attention on the so-called small- 
world phenomenon and a model proposed by Newman and 
Strogatz to replicate this phenomenon. 

Once the network is modeled, one is usually interested in 
two types of problems. The first involves structural prop- 
erties of the model. The second involves the performance 
of dynamical processes run on those networks. In the latter 
direction, the performance of random walks [12], Markov 
processes [6], gossip algorithms [5], consensus in a network 
of agents [16], [10], or synchronization of oscillators [21], 
[17], are very well reported in the literature. These dynamical 
processes are mostly studied in the traditional context of 
deterministic networks of relatively small size and/or regular 
structure. Even though many noteworthy results have been 
achieved for large-scale probabilistic networks [13]- [2], 
there is substantial reliance on numerical simulations. 

The eigenvalue spectrum of an undirected graph contains 
a great deal of information about structural and dynamical 
properties [7]. In particular, we focus our attention on the 
spectrum of the (combinatorial) Laplacian matrix uniquely 
associated with an undirected graph [3]. This spectrum 
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contains useful information about, for example, the number 
of spanning trees, or the stability of synchronization of a 
network of oscillators. We analyze the low-order moments 
of the Kirchhoff matrix spectrum corresponding to small- 
world networks. 

The paper is organized as follows. In Section II, we 
review the master stability function approach. In Section 
III, we derive closed-form expressions for the low-order 
moments of the Laplacian eigenvalue distribution associated 
with a probabilistic small-world network. Our expressions 
are valid for networks of asymptotically large size. Section 
IV applies our results to the problem of synchronization 
of a probabilistic small-world network of oscillators. The 
numerical results in this section corroborate our predictions. 

II. Synchronization of Nonlinear Oscillators 

In this section we review the master-stability-function 
(MSF) approach, proposed by Pecora and Carrol in [17], 
to study local stability of synchronization in networks of 
nonlinear oscillators. Using this approach, we reduce the 
problem of studying local stability of synchronization to the 
algebraic problem of studying the spectral support of the 
Laplacian matrix of the network. First, we introduce some 
needed graph-theoretical background. 

A. Spectral Graph Theory Background 

In the case of a network with symmetrical connections, 
undirected graphs provide a proper description of the net- 
work topology. An undirected graph G consists of a set of N 
nodes or vertices, denoted by V = {vi,...,v„}, and a set of 
edges E, where E 6 V x V. In our case, (v,-,Vy) 6 E implies 
(vj,Vj) €E, and this pair corresponds to a single edge with no 
direction; the vertices v, and v, are called adjacent vertices 
(denoted by v,- ~ vj) and are incident to the edge (v/,v/). We 
only consider simple graphs (i.e., undirected graphs that have 
no self-loops, so v, ^ vj for an edge (v,-,v ; ), and no more 
than one edge between any two different vertices). A walk 
on G of length k from vo to is an ordered set of vertices 
(vo,vi,...,vjt) such that (v;,V; + i) e E, for i = 0, 1, ...,k— 1; if 
v k = Vo the walk is said to be closed. 

The degree dt of a vertex v, is the number of edges incident 
to it. The degree sequence of G is the list of degrees, usually 
given in non-increasing order. The clustering coefficient, 
introduced in [22], is a measure of the number of triangles in 
a given graph, where a triangle is defined by the set of edges 
{{hi) -i U>k) ' suc h m at i <~ j <~ k ~ /. Specifically, we 

define clustering as the total number of triangles in a graph, 



T (G) , divided by the number of triangles in a complete (all- 
to-all) graph with N vertices, i.e., the coefficient is equal to 
T(G)/(»). 

It is often convenient to represent graphs via matrices. 
There are several choices for such a representation. For 
example, the adjacency matrix of an undirected graph G, 
denoted by A(G) = [a,- ; -], is defined entry-wise by ajj = 1 
if nodes ; and j are adjacent, and au = otherwise. (Note 
that an = for simple graphs.) Notice also that the degree 
dj can be written as dj = Ylj=i a ij- We can arrange the 
degrees on the diagonal of a diagonal matrix to yield the 
degree matrix, D = diag{dj). The Laplacian matrix (also 
called Kirchhoff matrix, or combinatorial Laplacian matrix) 
is defined in terms of the degree and adjacency matrices 
as L(G) = D{G) —A(G). For undirected graphs, L{G) is a 
symmetric positive semidefinite matrix [3]. Consequently, it 
has a full set of N real and orthogonal eigenvectors with real 
non-negative eigenvalues. Since all rows of L sum to zero, it 
always admits a trivial eigenvalue X\ = 0, with corresponding 
eigenvector vj = (1, 1, l) T ■ 

The moments of the Laplacian eigenvalue spectrum are 
central to our paper. Denote the eigenvalues of our N x N 
symmetric Laplacian matrix L(G) by = Ai (G) < ... < 
Ajv(G). The empirical spectral density (ESD) of L{G) is 
defined as 

MG (A) = if 5(A-A,), 

where 5(-) is the Dirac delta function. The A:-th order moment 
of the ESD of L(G) is defined as: 

**(G) = ^2>(G)* 

(which is also called the k-th order spectral momenQ). 

In the following subsection, we illustrate how a network 
of identical nonlinear oscillators synchronizes whenever the 
Laplacian spectrum is contained in a certain region on 
the real line. This region of synchronization is exclusively 
defined by the dynamics of each isolated oscillator and the 
type of coupling [17], [11]. This simplifies the problem of 
synchronization to the problem of locating the Laplacian 
eigenvalue spectrum. 

B. Synchronization as a Spectral Graph Problem 

Several techniques have been proposed to analyze the 
synchronization of coupled identical oscillators. In [23], 
well-known results in control theory, such as the passivity 
criterion, the circle criterion, and a result on observer design 
are used to derive synchronization criteria for an array 
of identical nonlinear systems. In [19], the authors use 
contraction theory to derive sufficient conditions for global 
synchronization in a network of nonlinear oscillators. We 
pay special attention to the master-stability-function (MSF) 

'Given that our interest is in networks of growing size (i.e., number of 
nodes N), a more explicit notation for p. and would perhaps have been 
fl^ and qi . However, for notational simplicity, we shall omit reference 
to N in there and other quantities in this paper. 



approach, [17]. This approach provides us with a criterion 
for local stability of synchronization based on the numerical 
computation of Lyapunov exponents. Even though quite 
different in nature, the mentioned techniques emphasize the 
key role played by the graph eigenvalue spectrum. 

In this paper we consider a time-invariant network of 
N identical oscillators, one located at each node, linked 
with 'diffusive' coupling. The state equations modeling the 
dynamics of the network are 

N 

x,=f(x,-) + rL fl '7 r ( x /- x ')> i = l,...,N (1) 

where x, represents an n-dimensional state vector corre- 
sponding to the 2-th oscillator. The nonlinear function f (•) 
describes the (identical) dynamics of the isolated nodes. The 
positive scalar y can be interpreted as a global coupling 
strength parameter. The « x « matrix T represents how 
states in neighboring oscillators couple linearly, and atj are 
the entries of the adjacency matrix. By simple algebraic 
manipulations, one can write down Eq. (Q~|i in terms of the 
Laplacian entries, L(G) = [/,-;], as 

N 

ii = f (x,) - y£ I.jTxj, for i = l,...,N. (2) 

7=1 

We say that the network of oscillators is at a synchronous 
equilibrium if X\(t) = X2(t) = ... = x N (t) = <j> (r ), where 
(f) represents a solution for x = f (x). In [17], the authors 
studied the local stability of the synchronous equilibrium. 
Specifically, they considered a sufficiently small perturbation, 
denoted by £;(?), from the synchronous equilibrium, i.e., 

Xi(f) = 0(0 + £;(')• 

After appropriate linearization, one can derive the following 
equations to approximately describe the evolution of the 
perturbations: 

e, = Df (/) ei(t) - y£ UjTejit), for i = l,...,N. (3) 

7=1 

where Df(f) is the Jacobian of f( ) evaluated along the 
trajectory (j) (t). This Jacobian is an n x « matrix with time- 
variant entries. Following the methodology introduced in 
[17], Eq. © can be similarity transformed into a set of linear 
time-variant (LTV) ODEs of the form: 

4= [Df (t) + (yA, (G)) T] for i = l,...,N, (4) 

where {A; (G)}i<,-<jv is the set of eigenvalues of L(G). Based 
on the stability analysis presented in [17], the network of 
oscillators in (Q]) presents a locally stable synchronous equi- 
librium if the corresponding maximal nontrivial Lyapunov 
exponents of (01 is negative for i = 2,...,N. 

Inspired in Eq. ©, Pecora and Carroll studied in [17] 
the stability of the following parametric LTV-ODE in the 
parameter a: 

|=[Df(0 + ar]S, (5) 



where Df (f) is the linear time-variant Jacobian in Eq. (|3). 
The master stability function (MSF), denoted by F(a), is 
defined as the value of the maximal nontrivial Lyapunov 
exponent of (0 as a function of a. Note that F (a) de- 
pends exclusively on f (•) and F, and is independent of the 
coupling topology, i.e., independent of L(G). The region of 
synchronization is, therefore, defined by the range of a > 
for which F (a) < 0. For a broad class of systems, the 
MSF is negative in the interval a 6 [0, <7 max ] = S (although 
more generic stability sets are also possible, we assume, 
for simplicity, this is the case in subsequent derivations). In 
order to achieve synchronization, the set of scaled nontrivial 
Laplacian eigenvalues, {Y^i}2<i<N> must be located inside 
the region of synchronization, S. This condition is equivalent 
to: yA 2 > and yk N < <7 max . 

We illustrate how to use of the above methodology in the 
following example: 
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Fig. 1. Periodic trajectory with period T = 5.749 in a Rossler oscillator 
when the parameters in Eqn. |6) take the values a = 0.2, b = 0.2, and c = 2.5. 



Example 1: Study the stability of synchronization of a 
ring of 6 coupled Rossler oscillators [14]. The dynamics of 
each oscillator is described by the following system of three 
nonlinear differential equations: 

M = -(yi + Zi), 
ji = Xi + ayt, 

Zi = b + Zi (xi - c) . 
The adjacency entries, a, 7 -, of a ring graph of six nodes 
are a,-; = 1 if j E 1) mod 6, (i— 1) mod 6}, for i = 
1,2,..., 6, and a, 7 - = otherwise. The dynamics of this ring 
of oscillators are defined by: 
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where we have chosen to connect the oscillators through 
their xi states exclusively. Our choice is reflected in the 
structure of the 3x3 matrix, Y, inside the summation in 
Eqn. ©. 

Numerical simulations of an isolated Rossler oscillator 
unveil the existence of a periodic trajectory with period 
T = 5.749 when the parameters in Eqn. ^ take the values 
a = 0.2, b = 0.2, and c = 2.5 (see Fig. 1). We denote this 
periodic trajectory by 0(f) = [(j) x (t) , (j) y (t) , (j) z (t)]. In our 
specific case, the LTP differential equation (0 takes the 
following form: 
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where the leftmost matrix in the above equation represents 
the Jacobian of the isolated Rossler evaluated along the 
periodic trajectory <j> (t) , and the rightmost matrix represents 

r. 

In Fig. 2, we plot the numerical values of the maximum 
Floquet exponent of Eqn. © for a € [0, 15], discretizing at 
intervals of length 0.2. This plot shows the range in which 
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Fig. 2. Numerical values of the maximum Floquet exponent of Eqn. (7) 
for o" £ [0, 15], discretizing at intervals of length 0.2. 



the maximal Floquet exponent is negative. This range of 
stability is S = (0,(7*), for a* re 4.7. The MSF criterion 
introduced in [17] states that the synchronous equilibrium 

is locally stable if the set of values {y A, (G)}, = 2 „ lies 

inside the stability range, S. For the case of a 6-ring 
configuration, the eigenvalues of L(G) are {0,1,1,3,3,4}, 

so the set {yA,},=2 „ is {y,y,3y,3y,4y}. Therefore, we 

achieve stability for yG (0, <7*/A„ (G)), where in our case 
a*/l n (G)^ 1.175. 

We now illustrate this result with several numerical sim- 
ulations. First, we plot in Fig. 3. a the temporal evolution of 
the Xi states of the 6-ring when y = 1.0. Observe how, since 
y€ (0,1.175), we achieve asymptotic synchronization. On 
the other hand, if we choose y= 1.3 ^ (0,1.175), the time 
evolution of the set of oscillators does not converge to a 
common trajectory (see Fig. 3.b); instead, the even and odd 
nodes settle into two different trajectories. 

In the next subsection, we propose an approach to estimat- 
ing the support of the eigenvalue distribution of large-scale 
probabilistic networks from low-order spectral moments. 
This allows us to predict synchronization in a large-scale 
Chung-Lu network. 
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Fig. 3. In Fig. a, we plot the temporal evolution of the X\ states of the 6-ring when 7= 1.0. In Fig.b, we plot the time evolution of the set of oscillators 
for 7 = 1.3 i (0,1.175). 



III. Spectral Analysis of Small- World 
Networks 

In this section we study the Laplacian eigenvalue spectrum 
of a variant of Watts -Strogatz small-world network [22]. 
After describing the model, we use algebraic graph theory to 
compute explicit expressions for the Laplacian moments of 
a small-world network as a function of its parameters. Our 
derivations are based on a probabilistic analysis of the ex- 
pected spectral moments of the Laplacian for asymptotically 
large small-world networks. 

A. Small-World Probabilistic Model 

We consider a one-dimensional lattice of N vertices, 
{v\,...,vn}, with periodic boundary conditions, i.e., on 
a ring, and connect each vertex v,- to its 2k clos- 
est neighbors, i.e., V; is connected to the set of nodes 
{vj-.je[(i-k)modN,(i + k)modN]}. Then, instead of 
rewiring a fraction of the edges in the regular lattice as 
proposed by Watts and Strogatz [22], we add some random 
'shortcuts' to the one-dimensional lattice. These shortcuts are 
added by independently assigning edges between each pair of 
nodes , 1 < i < j <N with probability p. The resulting 
small-world graph is intermediate between a regular lattice 
(achieved for p = 0) and a classical random graph (achieved 
for p = 1). In general, small-world networks share properties 
with both the regular grid and the classical random graph for 
< p < 1 . In particular, they show the following apparently 
contradictory features: 

(i) most nodes are not neighbors of one another (such as 
in a regular grid), and 

(ii) most nodes can be reached from every other node by 
a small number of steps (such as in a random graph). 

An interesting property observed in this model was the 
following: for small probability of rewiring, p< 1, the 
number of triangles in the network is nearly the same as that 
of the regular lattice, but the average shortest-path length is 
close to that of classical random graphs. In the rest of the 
paper we assume we are in the range of p in which this 
property holds, in particular, we will prescribe p to be r/N, 
for a given parameter r. 

In the coming sections we shall study spectral properties 
of the Laplacian matrix associated to the above small-world 



model. In our derivations we will need the probabilistic 
distribution for the degrees. It is well known that, for asymp- 
totically large graphs, the degree distribution of a classical 
random graph with average degree r is a Poisson distribution 
with rate r. Hence, the degree distribution of the above small- 
world network is 



Pr(di = d) 



0, 

r d-2k e -r 

(d-2k) ! 



for d < 2k, 
for d > 2k, 



(8) 



which corresponds to a Poisson with parameter r 'shifted' 2k 
units. The Poisson distribution is shifted to take into account 
the degree of the regular 2£-neighbors ring superposed to the 
random shortcuts. 

Furthermore, it is well known that the clustering coeffi- 
cient (or, equivalently, the number of triangles) of the regular 
2A:-neighbors rings is very lightly perturbed by the addition 
of random shortcuts for p = r/N. In particular, one can prove 
the following result: 



1 /2k 
L[T] = (l+o(l))-N( 



where the dominant term, hN(^ 2 ), corresponds to the exact 
number of triangles in a 2£-neighbors ring with N nodes. 

In the following section, we shall derive explicit ex- 
pressions for the first low-order spectral moments of the 
Laplacian matrix associated with the small-world model 
herein described. Even though our analysis is far from 
complete, in that only low-order moments are provided, 
valuable information regarding spectral properties can be 
retrieved from our results. 

B. Algebraic Analysis of Spectral Moments 

In this section we deduce closed-form expressions for the 
first three moments of the Laplacian spectrum of any simple 
graph G. First, we express the spectral moments as a trace 
using the following identity: 



(9) 



1 N 1 



<7*(G) = -£A ( -(G)*=-trL(G) 



N 



(10) 



This identity is derived from the fact that trace is conserved 
under diagonalization (in general, under any similarity trans- 
formation). In the case of the first spectral moment, we obtain 



1 

N' 



-tr(D-A) 



N 

1=1 



where d is the average degree of the graph. For analytical 
and numerical reasons, we define the normalized Kirchhoff 
moment as 

2* = i£(W=4^tr(Z>-A)*. (11) 



The fact that D and A do not commute forecloses the pos- 
sibility of using Newton's binomial expansion on (D—A), 
On the other hand, the trace operator allows us to cyclically 
permute multiplicative chains of matrices. For example, 
tr (AAD) =tr (ADA) =tr (DAA). Thus, for words of length 
k < 3, one can cyclically arrange all binary words in the 
expansion of ( fTTT i into the standard binomial expression: 



Ik 



k 

E 

a=0 



^4^tr (A a D k - 



for k < 3. (12) 



Also, we can make use of the identity tr(A a D 
E*Li(A a ) fl </?~ a to write 



,k-a\ 



k N 
a=0i=l 



( 



1 \ a 



d k N 



for it < 3. (13) 



Note that this expression is not valid for k > 4. For example, 
for k — 4, we have that tr (AADD) ^tv(DADA) . 

We now analyze each summand in expression ( [TBI from 
a graph-theoretical point of view. Specifically, we find a 
closed-form solution for each term tr(A'D ; ), for all pairs 
1 < i + j < 3, as a function of the degree sequence and the 
number of triangles in the network. In our analysis, we make 
use of the following results from [3]: 

Lemma 2: The number of closed walks of length a in a 
graph G, joining node i to itself, is given by the ;-th diagonal 
entry of the matrix A a . 

Corollary 3: Let G be a simple graph. Denote by f, the 
number of triangles touching node ;. Then, 



(A)„. = 0, (A ) .. = di, and (A 3 ), = 2f,-. 



(14) 



After substituting ( fT~4T > into (fT3l , and straightforward al- 
gebraic simplifications, we obtain the following exact ex- 
pression for the low-order normalized spectral moments of 
a given Kirchhoff matrix K: 



1k = 



1. 

l 

Ay 2 

Nd 3 



[(Xl^ + 311^-67] 



for k = l, 
for k = 2, 

for k= 3, 



where T = jI-L^; is the total number of triangle^ in the 
network. 

It is worth noting how our spectral results are written in 
terms of two widely reported measurements, [15]: the degree 
sequence and the clustering coefficient (which provides us 
with the total number of triangles.) This allows us to compute 
low-order spectral moments of many real-world networks 
without performing an explicit eigenvalue decomposition. 



C. Probabilistic Analysis of Spectral Moments 

In this section, we use Eq. ( fTBT ) to compute the first 
three expected Laplacian moments of the small-world model 
under consideration. The expected moments can be computed 
if we had explicit expressions for the moments of the 
degrees, E[rf,], E[c/?], and E[rf?], and the expected number 
of triangles, K[T]. Since we know the degree distribution (H)) 
for this model, the moments of the degrees can be computed 
to be: 



E[di] 
E[df] 
E[df] 



r + 2k 

2 



(16) 



r z + (l+4k)r + 4-k 2 , 
r 3 + l5{3 + 6k)r 2 +(l+6k+\2k 2 ) r-8k\ 



We can therefore substitute the expressions (|9]l and ( TToT l in 
Eq. $15[ in order to derive the following expressions for the 
(non-normalized) expected Laplacian moments for N — > °°: 



E[qi] 
E[?2] 

E[q 3 ] 



r + 2k, 

r 2 + (4k + 2)r + 4k 2 



(17) 



■2k, 



r 3 + (6k + 6) r 2 + ( 12k 2 + 18Jt + 4) r 



-M 2 + U i + 2k. 



In the following table we compare the numerical values 
of the Laplacian moments corresponding to one random 
realization of the model under consideration with the an- 
alytical predictions in (T% . In particular, we compute the 
moments for a network of = 512 nodes with parameters 
p = r/N =4/N and k = 3. It is important to point out that the 
indicated numerical values are obtained for one realization 
only, with no benefit from averaging. 



Moment order 


I st 


2 nd 


yd 


Numerical realization 


10.14 


116.96 


1,467.6 


Analytical expectations 


10 


114 


1,431 


Relative error 


1.38% 


2.53% 


2.49% 



(15) 



In the next subsection, we use an approach introduced in 
[18] to estimate the support of the eigenvalue distribution 
using the first three spectral moments. In coming sections, 
we shall use this technique to predict whether the Laplacian 
spectrum lies in the region of synchronization. 

2 A triangle is defined by a set of (undirected) edges {(i,j),(j,k),(k,i)} 
such that i ~ j ~ k ~ i. 



D. Piecewise-Linear Reconstruction of the Laplacian Spec- 
trum 

Our approach, described more fully in [18], approximates 
the spectral distribution with a triangular function that ex- 
actly preserves the first three moments. We define a triangular 
distribution t (A) based on a set of abscissae x\ <x 2 <x 3 as 



t(A):= 



i£r(A-i). 

0, 



for A e [x\,x 2 ), 
for A € [x 2 ,x 3 ] , 
otherwise. 



where h =2/{x^—X\). The first three moments of this 
distribution, as a function of the abscissae, are given by 

1 



M 2 
M 3 



(18) 



£ 3 +x\X2 -I-JC1JC3 +X2X3J 



— (xf +x\x2 +x}x3 +x\ +x\x\ 



■ X2X3 ■ 



-4- 



2 2 \ 

-JC3JC1 +X3X2 +X\X2Xt> J 



Our task is to find the set of values {xi ,x 2 ,X3} in order to fit 
a given set of moments {M\,M2,Mj,}. The resulting system 
of algebraic equations is amenable to analysis, based on the 
observation that the moments are symmetric polynomial^. 
Following the methodology in [18], we can find the abscissae 
{xi,X2,x$} as roots of the polynomial: 



where 



^-nix 2 +n 2 x-n 3 = o, (19) 



IIi = 3 Mi, (20) 

n 2 = 9M?-6M 2 , 

n 3 = 27Mj , -36MiM 2 + 10M 3 . 

The following example illustrates how this technique pro- 
vides a reasonable estimation of the Laplacian spectrum for 
small-world Networks. 

Example 4: Estimate the spectral support of the small- 
world model described in Subsection IIII-AI for parameters 
N = 512, p = 4/N and k = 3. In subsection IlITCl we 
computed the expected spectral moments of this particular 
network to be {Mi = 10, M 2 = 114,M 3 = 1,431}. Thus, we 
apply the above technique with these particular values of the 
moments to compute the following set of abscissae for the tri- 
angular reconstruction {x\ = 1.577, x 2 = 8.662,X3 = 19.76}. 
In Fig. 4 we compare the triangular function that fits the 
expected spectral moments with the histogram of the eigen- 
values of one random realization of the Laplacian matrix. We 
also observe that any random realization of the eigenvalue 
histograms of the Laplacian is remarkably close to each 
other. Although a complete proof of this phenomenon is 
beyond the scope of this paper, one can easily proof using 
the law of large numbers that the distribution of spectral 
moments in < TT~5T > concentrate around their mean values. 

3 A symmetric polynomial on variables (xi,X2,X3) is a polynomial that is 
unchanged under any permutation of its variables. 




Fig. 4. Comparison between the histogram of the eigenvalues of one 
random realization of the Laplacian matrix of a small-world model with 
parameters N = 512, p = 4/N and k = 3, and the triangular function that 
fits the expected spectral moments. 




Fig. 5. Comparison between the values of the triangular abscissae X\ and 
X3 with the extreme points of the Laplacian spectral support, A2 and A„, for 
a small-world network with N = 512 nodes, k = 3, and p in the range of 
values [1/JV : 0.01/JV : 10/JV] . 



We observe that the above estimation is valid for a large 
range in the values of the parameters. For example, in Fig. 
5, we compare the values of the triangular abscissae x\ and 
x 3 with the extreme points of the Laplacian spectral support, 
A 2 and A„, for a small-world network with N = 512 nodes, 
k = 3, and p in the range of values [1 /N : 0.01 /N : 10/N] . It 
is important to point out that, in this case too, the numerical 
values for the eigenvalues are obtained for one realization 
only, with no benefit from averaging. In the next section, 
we propose a methodology which uses results presented 
in previous sections to predict the local stability of the 
synchronous state in a small-world network of oscillators. 

IV. Analytical Estimation of Synchronization 

In this section we use the expressions in ( TTTb and the 
triangular reconstruction in the above subsection to predict 
synchronization in a large small-world network of coupled 
nonlinear oscillators. Specifically, we study a network of 
coupled Rossler oscillators, as those in Example Q] We build 
our prediction based on the following steps: 

1) Determine the region of synchronization following 
the technique presented in Subsection III-BI As illus- 
trated in Example Q] the set of scaled eigenvalues 
{yX i }i=2 n must lie in a certain region of stability, 



S, to achieve synchronization (in our example S = 
(0,ff*w4.7)). 

2) Compute the expected spectral moments of the Lapla- 
cian eigenvalue spectrum for a given set of parameters 
using the set of Eqns. in (T7\ . 

3) Estimate the support of the Laplacian eigenvalue spec- 
trum, {A/ }i=2 ...jv> using the methodology presented 
in Subsection IIII-DI From Example 2J we have that 
s/ = 1.57 and s u = 19.76 are good estimates of the 
lower and upper extremes of the spectral support, 
respectively. 

4) Compare the region of stability in Step 1 with the 
estimation of the spectral support in Step 3, i.e., 
(1.57 7,19.76 y). 

Following the above steps, one can easily verify that our 
estimated spectral support, (1.57 7,19.76 7), lies inside the 
region of stability, (0,(7* «4.7), for < 7 < 4.7/19.76 w 
0.238. Therefore, the small-world network of 512 coupled 
Rossler oscillators is predicted to synchronize whenever the 
global coupling strength satisfies 76 (0,0.238). 

A. Numerical Results 

In this section we present numerical simulations support- 
ing our conclusions. We consider a set of identical 512 
Rossler oscillators (as the one described in Example Q]) 
interconnected through the Small-World network defined 
in Example (p = 4/N and k = 3). Using the methodology 
proposed above, we have predicted that the synchronous state 
of this system is locally stable if the coupling parameter 7 
lies in the interval (0,0.238). We run several simulations with 
the dynamics of the oscillators presenting different values of 
the global coupling strength 7. For each coupling strength, 
we present two plots: (i) the evolution of the 5 12 x-states of 
the Rossler oscillators in the time interval < t < 40, and ( ii) 
the evolution of Xi (t)—x(t) for all i, where x(t) = ^ Y,i x i (0- 
Since our stability results are local, we have to carefully 
choose the initial states for the network of oscillators. For 
our particular choice of parameters, the (isolated) Rossler 
oscillator presents a stable limit cycle. For our simulations, 
we have chosen as initial condition for each oscillator in 
the network a randomly perturbed version of a particular 
point of this stable limit cycle. This particular point is So = 
(3.5119,-3.5332,0.2006). We have chosen the perturbed 
initial state for the i-th oscillator to be So + e,, where e, is a 
uniformly distributed random variable in the 3 -dimensional 
cube [— 2,2] 3 , and e,- is independent of e ; for i ^ j.l 

In our first simulation, we use a coupling strength 7 = 
0.1 G (0,0.238); thus, we predict the synchronous state to 
be locally stable. Fig. 6 (a) and (b) represents the dynamics 
A-states for the 512 oscillators in the small-world network. 
In this case, we observe a clear exponential convergence of 
the errors to zero. In the second simulation, we choose 7 = 
0.3 ^ (0,0.238); thus, we predict the synchronous state to 
be unstable. In fact, we observe in Figs. 7. a and 7.b how 
synchronization is clearly not achieved. 



V. Conclusions and Future Research 

In this paper, we have studied the eigenvalue distribu- 
tion of the Laplacian matrix of a large-scale small-world 
networks. We have focused our attention on the low-order 
moments of the spectral distribution. We have derived ex- 
plicit expressions of these moments as functions of the 
parameters in the small-world model. We have then applied 
our results to the problem of synchronization of a network 
of nonlinear oscillators. Using our expressions, we have 
studied the local stability of the synchronous state in a large- 
scale small-world network of oscillators. Our approach is 
based on performing a triangular reconstruction matching the 
first three moments of the unknown spectral measure. Our 
numerical results match our predictions with high accuracy. 
Several questions remain open. The most obvious extension 
would be to derive expressions for higher-order moments of 
the Kirchhoff spectrum. A more detailed reconstruction of 
the spectral measure can be done based on more moments. 
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